This notebook is a fork of "AnkleSLIP - find minimal model". Here, only a the "SLIP parameter prediction and model self-consistency check" figure is created last edits:
In [1]:
# run this to connect an ipython qtconsole to the notebook
%connect_info
In [1]:
%cd
%cd mmnotebooks
In [2]:
# system libs
import os
import sys
import re
from copy import deepcopy
# shai's lib
import libshai.util as ut
# my libs
import mutils.io as mio
import mutils.misc as mi
import mutils.statistics as st
import mutils.FDatAn as fda
from mutils.io import build_dataset
import models.sliputil as su
import models.slip as sl
from models.sliputil import getControlMaps
from models.sliputil import get_auto_sys
import mutils.plotting as mplt
conf = mio.saveable()
# define and import data
conf.subject = 2 #available subjects: 1,2,3,7. Other subjects have comparatively bad data quality
conf.ttype = 1 # 1: free running, 2: metronome running (data incomplete!)
conf.n_factors = 5 # 1-5 factors to predict SLIP parameters
conf.n_factors_doc = "how many (optimal) factors to select of the full kinematic state"
conf.exclude_IC_from_factors = False # exclude the initial conditions from kinematic state
# detrending options
conf.dt_window = 30
conf.dt_medfilter = False # use median filter instead of mean
# select ankle-SLIP instead of "factors without CoM state" ?
conf.select_ankle_SLIP = True
conf.cslip_forceZeroRef = True # reference values for controlled SLIP maps must be zero or not
conf.cslip_forceZeroRef_doc = "reference values for controlled SLIP maps must be zero or not"
# to compute periodic SLIP orbit: average over parameters or initial conditions?
conf.po_average_over_IC = True
conf.po_average_over_IC_doc = "average over IC's and T,ymin (alt: parameters) for reference SLIP"
# startup
conf.startup_compute_full_maps = True # compute return maps on data loading?
conf.startup_n_full_maps = 40 # how many full-stride maps to compute for averaging
conf.startup_compute_PCA = False # compute PCA on startup?
conf.display()
In [3]:
%%writefile tmp/step1.py
# load kinematic data
ws1 = mio.saveable()
#ws1.k = mio.KinData(data_dir = '/home/moritz/data/2011-mmcl_mat/')
ws1.k = mio.KinData(data_dir = 'data/2011-mmcl_mat/')
ws1.k_doc = "generic kinematic data access object (empty selection)"
ws1.k.load(conf.subject, conf.ttype)
ws1.k_doc = "KinData object: s:" + str(conf.subject) + ", t:" + str(conf.ttype) + ", selection: ankles"
ws1.k.selection = ['r_anl_y - com_y', 'l_anl_y - com_y', 'com_z']
ws1.SlipData = [mi.Struct(mio.mload('data/2011-mmcl_mat/SLIP/new/params3D_s%it%ir%i.dict' %
(conf.subject, conf.ttype, rep)))
for rep in ws1.k.reps]
ws1.SlipData_doc = "SlipData for s:" + str(conf.subject) + ", t:" + str(conf.ttype)
ws1.display()
# huge data selection, but not every dimension of every marker
ws1.full_markers = [ 'com_x', 'com_y', 'com_z', 'r_mtv_x - com_x',
'r_anl_x - com_x', 'r_kne_x - com_x',
'r_sia_x - com_x',
'l_mtv_x - com_x',
'l_anl_x - com_x', 'l_kne_x - com_x', 'l_sia_x - com_x',
'sacr_x - com_x',
'cvii_x - com_x', 'r_wrl_x - com_x',
'r_elb_x - com_x', 'l_elb_x - com_x',
'l_wrl_x - com_x',
'r_anl_y - com_y', 'r_kne_y - com_y', 'r_trc_y - com_y',
'r_sia_y - com_y',
'sacr_y - com_y', 'l_anl_y - com_y',
'l_kne_y - com_y', 'l_trc_y - com_y',
'l_sia_y - com_y',
'cvii_y - com_y',
'r_hea_x - com_x', 'l_hea_x - com_x', 'r_hea_y - com_y', 'l_hea_y - com_y',
'r_hea_z - com_z', 'l_hea_z - com_z',
'r_elb_y - com_y',
'r_acr_y - com_y', 'l_acr_y - com_y',
'l_elb_y - com_y', 'r_mtv_z - com_z',
'r_anl_z - com_z',
'r_trc_z - com_z',
'l_mtv_z - com_z',
'l_anl_z - com_z', 'l_trc_z - com_z',
'sacr_z - com_z',
'r_wrl_z - com_z',
'r_acr_z - com_z', 'l_acr_z - com_z',
'l_wrl_z - com_z']
ws1.full_markers_doc = 'representative selection of markers for quasi-full state space information'
print "len of 'ws1.full_markers':", len(ws1.full_markers), " (all: 84)"
ws1.k.selection = ws1.full_markers #all_markers
ws1.dataset_full = build_dataset(ws1.k, ws1.SlipData, dt_window=conf.dt_window, dt_median=conf.dt_medfilter)
ws1.dataset_full_doc = ("reference 'dataset' for s:" + str(ws1.k.subject) + " t:" + str(ws1.k.ttype) +
" with 'full' selection")
# "normalized" dataset: velocities scaled by small factor
ws1.dataset_full.n_kin_r = ws1.dataset_full.all_kin_r.copy()
ws1.dataset_full.n_kin_l = ws1.dataset_full.all_kin_l.copy()
ws1.dataset_full.n_kin_r[:, len(ws1.k.selection)-2:] /= 11.
ws1.dataset_full.n_kin_l[:, len(ws1.k.selection)-2:] /= 11.
ws1.dataset_full.n_kin_r -= mean(ws1.dataset_full.n_kin_r, axis=0)
ws1.dataset_full.n_kin_l -= mean(ws1.dataset_full.n_kin_l, axis=0)
ws1.dataset_full.n_kin_r = fda.dt_movingavg(ws1.dataset_full.n_kin_r, conf.dt_window, conf.dt_medfilter)
ws1.dataset_full.n_kin_l = fda.dt_movingavg(ws1.dataset_full.n_kin_l, conf.dt_window, conf.dt_medfilter)
ws1.dataset_full.n_kin_r_doc = "detrended 'all_kin_r'; velocity has been scaled by factor 1/11"
ws1.dataset_full.n_kin_l_doc = "detrended 'all_kin_l'; velocity has been scaled by factor 1/11"
In [4]:
%%writefile tmp/step2.py
## prepare data
# *NOTE* the first three labels have to be "com_x", "com_y", "com_z". If you edit this, make sure to edit
# the cell where data are further processed...
ws1.k.selection = ws1.full_markers
ws1.k_doc = "KinData object: s:" + str(ws1.k.subject) + ", t:" + str(ws1.k.ttype) + ", selection: 'full'"
ws1.dataset = ws1.dataset_full
ws1.dataset_doc = ws1.dataset_full_doc
# another consistency check
stepnr = 698 # select any stride between 0 and ~1800 (subject dependend)
mass = mean([d.mass for d in ws1.SlipData])
# su.finalState actually performs a one-step integration of SLIP
print "consistency check:"
print "------------------"
print "simres [right step]:", su.finalState(ws1.dataset.all_IC_r[stepnr, :], ws1.dataset.all_param_r[stepnr, :],
{'m' : mean(ws1.dataset.masses), 'g' : -9.81})
print "subsequent left apex:", ws1.dataset.all_IC_l[stepnr, :]
print "==========="
print "simres [left step]:", su.finalState(ws1.dataset.all_IC_l[stepnr, :], ws1.dataset.all_param_l[stepnr, :],
{'m' : mean(ws1.dataset.masses), 'g' : -9.81})
print "subsequent right apex:", ws1.dataset.all_IC_r[stepnr + 1, :]
In [5]:
%%writefile -a tmp/step2.py
# Here, the main predictors ("directions") are computed. This is computationally quite fast
if conf.exclude_IC_from_factors:
# indices that do *not* represent CoM (only relative positions, indicated by "minus" sign)
sel_idx = array([nr for nr, label in enumerate(ws1.dataset.kin_labels) if '-' in label])
ws1.facs_r_doc = ("directions of 'factors' for a right step excluding CoM, s:" +
str(ws1.k.subject) + ", t:" + str(ws1.k.ttype))
ws1.facs_l_doc = ("directions of 'factors' for a left step excluding CoM, s:" +
str(ws1.k.subject) + ", t:" + str(ws1.k.ttype))
else:
sel_idx = arange(len(ws1.dataset.kin_labels))
ws1.facs_r_doc = ("directions of 'factors' for a right step including CoM, s:" +
str(ws1.k.subject) + ", t:" + str(ws1.k.ttype))
ws1.facs_l_doc = ("directions of 'factors' for a left step including CoM, s:" +
str(ws1.k.subject) + ", t:" + str(ws1.k.ttype))
ws1.facs_labels = [ws1.dataset.kin_labels[x] for x in sel_idx]
ws1.facs_labels_doc = "kinematic labels for the components of the 'factors'"
ws1.facs_r = st.find_factors(ws1.dataset.s_kin_r[:, sel_idx].T, ws1.dataset.s_param_r.T, k=conf.n_factors)
ws1.facs_l = st.find_factors(ws1.dataset.s_kin_l[:, sel_idx].T, ws1.dataset.s_param_l.T, k=conf.n_factors)
ws1.fscore_r = dot(ws1.facs_r.T, ws1.dataset.s_kin_r[:, sel_idx].T).T
ws1.fscore_r_doc = "scores of data when projected on the factors (right)"
ws1.fscore_l = dot(ws1.facs_l.T, ws1.dataset.s_kin_l[:, sel_idx].T).T
ws1.fscore_l_doc = "scores of data when projected on the factors (left)"
ws1.raw_factor_space_r = ws1.dataset.s_kin_r[:, sel_idx]
ws1.raw_factor_space_l = ws1.dataset.s_kin_l[:, sel_idx]
ws1.raw_factor_space_l_doc = "part of state space used for computing the 'factors'"
ws1.raw_factor_space_l_doc = "part of state space used for computing the 'factors'"
# ---------------------visualize
figure(figsize=(20,3))
pcolor(ws1.facs_r.T)
gca().set_xticks(arange(len(ws1.facs_labels)) + .5)
gca().set_xticklabels(ws1.facs_labels, rotation=90)
grid('on')
title('magnitude of the factors')
gca().set_yticks(arange(ws1.facs_r.shape[1]) + .5)
gca().set_yticklabels(arange(ws1.facs_r.shape[1]) + 1)
gca().set_ylabel('factor #')
pass
In [6]:
%%writefile tmp/step3.py
# --- --- compute
def predTestFacs(ICr, ICl, fspacer, fspacel, pr, pl, nboot = 50):
"""
performs prediction test, computing the factors and keeping their direction constant in prediction
"""
# --- prepare output
rvpr = [] # relative variance for parameters (right step -> left(!) parameters)
rvfr = [] # relative variance for factor state (right step -> left(!) state)
rvffr = [] # relative variance for factor state (right step -> left(!) state)
# --- bootstrap loop
idat_fr = hstack([ICr, fspacer])
idat_fl = hstack([ICr, fspacel])
print "=" * nboot
for rep in range(nboot):
sys.stdout.write('.')
# split into regression and prediction (indices)
ridx = randint(0, idat_fr.shape[0], idat_fr.shape[0])
pidx = fda.otheridx(ridx, idat_fr.shape[0])
# compute factor direction
facsr = st.find_factors(fspacer[ridx, :].T, pr[ridx, :].T, k=conf.n_factors)
idat_facr = vstack([ICr.T, dot(facsr.T, fspacer.T)]).T
facsl = st.find_factors(fspacel[ridx, :].T, pl[ridx, :].T, k=conf.n_factors)
idat_facl = vstack([ICl.T, dot(facsl.T, fspacel.T)]).T
# compute maps
idatRr = idat_facr[ridx, :]
idatRfr = idat_fr[ridx, :]
odatRpr = pr[ridx, :]
odatRfr = idat_facl[ridx, :]
Ap_r = dot(odatRpr.T, pinv(idatRr.T))
Af_r = dot(odatRfr.T, pinv(idatRr.T))
Aff_r = dot(odatRfr.T, pinv(idatRfr.T))
# compute predictions
idatPr = idat_facr[pidx, :]
idatPfr = idat_fr[pidx, :]
odatPpr = pr[pidx, :]
odatPfr = idat_facl[pidx, :]
predpr = dot(Ap_r, idatPr.T).T
predfr = dot(Af_r, idatPr.T).T
predffr = dot(Aff_r, idatPfr.T).T
# compute relative remaining variance
rvpr.append(var(odatPpr - predpr, axis=0) / var(odatPpr, axis=0))
rvfr.append(var(odatPfr - predfr, axis=0) / var(odatPfr, axis=0))
rvffr.append(var(odatPfr - predffr, axis=0) / var(odatPfr, axis=0))
print "done"
return vstack(rvpr), vstack(rvfr), vstack(rvffr)
ICr, ICl, RFSr, RFSl, Pr, Pl = [fda.dt_movingavg(x, conf.dt_window, conf.dt_medfilter) for x in
[ws1.dataset.all_IC_r, ws1.dataset.all_IC_l, ws1.raw_factor_space_r, ws1.raw_factor_space_l,
ws1.dataset.s_param_r, ws1.dataset.s_param_l]]
# SLIP param prediction using factors, factors using factors, factors using FS
rvpr_fac, rvfr_fac, rvfr_full = predTestFacs(ICr, ICl, RFSr, RFSl, Pr, Pl, nboot=100)
rvpr_aug = st.predTest(hstack([ICr, roll(Pl, 1, axis=0)])[1:, :], Pr[1:, :], nboot=100)
rvpr_full = st.predTest(hstack([ICr, RFSr]), Pr, nboot=100)
raidx = array([idx for idx, lbl in enumerate(ws1.dataset_full.kin_labels) if 'r_anl' in lbl])
laidx = array([idx for idx, lbl in enumerate(ws1.dataset_full.kin_labels) if 'l_anl' in lbl])
rank = fda.dt_movingavg(ws1.dataset_full.all_kin_r[:, raidx], conf.dt_window, conf.dt_medfilter)
lank = fda.dt_movingavg(ws1.dataset_full.all_kin_l[:, laidx], conf.dt_window, conf.dt_medfilter)
rvpr_ank = st.predTest(hstack([ICr, rank]), Pr, nboot=100)
rvankr_ank = st.predTest(hstack([ICr, rank]), hstack([ICl, lank]) , nboot=100)
rvankr_full = st.predTest(hstack([ICr, RFSr]), hstack([ICl, lank]) , nboot=100)
print "done"
In [7]:
# --- compute predictions for all 4 subjects:
a_rvpr_full, a_rvpr_fac, a_rvpr_ank, a_rvpr_aug = [], [], [], []
a_rvfr_fac, a_rvfr_full, a_rvankr_ank, a_rvankr_full = [], [], [], []
for conf.subject in [1,2,3,7]:
print "=" * 20, "subject", conf.subject, "=" * 20
%run -i tmp/step1.py
%run -i tmp/step2.py
%run -i tmp/step3.py
# -> store computed variance reductions
a_rvpr_full.append(rvpr_full)
a_rvpr_fac.append(rvpr_fac)
a_rvpr_ank.append(rvpr_ank)
a_rvpr_aug.append(rvpr_aug)
a_rvfr_fac.append(rvfr_fac)
a_rvfr_full.append(rvfr_full)
a_rvankr_ank.append(rvankr_ank)
a_rvankr_full.append(rvankr_full)
In [9]:
rvfr_fac.shape
Out[9]:
In [11]:
a_rvpr_full[0].shape
Out[11]:
In [19]:
# --- visualize
# --- --- select subject for "full" visualization
vis_subj = 2 # 0,1,2,3
nvs = list(set(range(4)) - set([vis_subj]))
vis_subj_str = ['1', '2', '3', '7'][vis_subj]
rvpr_full = a_rvpr_full[vis_subj]
rvpr_fac = a_rvpr_fac[vis_subj]
rvpr_ank = a_rvpr_ank[vis_subj]
rvpr_aug = a_rvpr_aug[vis_subj]
rvfr_fac = a_rvfr_fac[vis_subj]
rvfr_full = a_rvfr_full[vis_subj]
rvankr_ank = a_rvankr_ank[vis_subj]
rvankr_full = a_rvankr_full[vis_subj]
# --- --- create figure
figure(figsize=(14,8))
# --- --- left subplot
ax0 = subplot(2,1,1)
posP = arange(rvpr_full.shape[1]) * 3
b1 = boxplot(rvpr_full, positions = posP + .1, widths=.30)
b2 = boxplot(rvpr_fac, positions = posP + .5, widths=.30)
b3 = boxplot(rvpr_ank, positions = posP + .9, widths=.30)
b4 = boxplot(rvpr_aug, positions = posP + 1.3, widths=.30)
mi.recolor(b1, '#006767')
mi.recolor(b2,'k')
mi.recolor(b3, '#009900')
mi.recolor(b4,'r')
fmt = {'marker' : '>',
'markersize' : 3,
'markeredgecolor' : None}
colors_ = mplt.colorset_distinct
for idx in nvs:
for x in arange(5):
plot(posP[x] + .1 - .1, median(a_rvpr_full[idx][:,x]), color=colors_[0], **fmt)
plot(posP[x] + .5 - .1, median(a_rvpr_fac[idx][:,x]), color=colors_[1], **fmt)
plot(posP[x] + .9 - .1, median(a_rvpr_ank[idx][:,x]), color=colors_[2], **fmt)
plot(posP[x] + 1.3 - .1, median(a_rvpr_aug[idx][:,x]), color=colors_[3], **fmt)
# --- --- right subplot
outliers = '' # enter '+' to show outliers
ax1 = subplot(2,1,2)
pos1 = hstack([arange(3) * 3, arange(conf.n_factors) * 1.8 + 10])
pos2 = hstack([arange(3) * 3, arange(6) * 1.8 + 18]) + .8
b1 = boxplot(rvfr_fac, sym=outliers, positions=pos1 + .1, widths=.20)
b2 = boxplot(rvfr_full, sym=outliers, positions=pos1 + .4, widths=.20)
ax1b = ax1.twiny()
b3 = boxplot(rvankr_ank, sym=outliers, positions=pos2 + .1, widths=.20)
b4 = boxplot(rvankr_full, sym=outliers, positions=pos2 + .4, widths=.20)
mi.recolor(b1,'k')
mi.recolor(b2,'k',2)
mi.recolor(b3,'#009900')
mi.recolor(b4,'#009900',2)
# remove outliers
ax0.set_xticks(arange(rvpr_full.shape[1]) * 3 + .7)
ax0.set_xticklabels(['$k$', r'$\alpha$', r'$\L_0$', r'$\beta$', r'$\Delta E$'])
ax0.set_xlabel('SLIP parameter')
ax0.set_ylabel('relative remaining variance')
ax0.plot([-1, rvpr_full.shape[1] * 3 - 1], [1, 1], 'k--')
ax0.set_xlim(-1, rvpr_full.shape[1] * 3 - 1)
ax0.set_ylim(0, 1.15)
xticks1 = pos1 + .45
xticks1[:3] += .25
ax1.set_xticks(xticks1)
xticks2 = pos2[3:] + .45
ax1b.set_xticks(xticks2)
ax1.grid('on')
ax1b.grid('on')
ax1.set_xticklabels(['com_z', 'v_com_x', 'v_com_y', 'f1', 'f2', 'f3', 'f4', 'f5'])
ax1b.set_xticklabels(['ank_x', 'ank_y', 'ank_z', 'v_ank_x', 'v_ank_y', 'v_ank_z'])
ax1.set_xlim(-1.5, pos2[-1] + 2)
ax1b.plot([-1.5, pos2[-1] + 2], [1, 1], 'k--')
ax1b.set_xlim(-1.5, pos2[-1] + 2)
ax1.set_ylim(0, 1.15)
ax1.set_ylabel('relative remaining variance')
ax1b.set_title('Reduced model inner consistency check. Boxplotes for subject {}\n\n'.format(vis_subj_str))
ax0.set_title('SLIP parameter prediction. Boxplots for subject {}'.format(vis_subj_str))
for idx in nvs:
for x in arange(len(pos1)):
plot([pos1[x] -.6, pos1[x] -.1],
[median(a_rvfr_fac[idx][:,x]), median(a_rvfr_full[idx][:,x])], '.-', color='k', linewidth=1)
for x in arange(len(pos2)):
plot([pos2[x] + .5, pos2[x] + 1],
[median(a_rvankr_ank[idx][:,x]), median(a_rvankr_full[idx][:,x])], '.-', color='g', linewidth=1)
subplots_adjust(hspace=.5)
#savefig('tmp/prediction_and_consistency.pdf')
savefig('tmp/prediction_and_consistency_f3.pdf')
pass #suppress output
In [9]:
# --- required imports and definitions
import matplotlib.font_manager as fmng
FM = fmng.FontManager()
fnt = FM.findfont('Arial')
fp = fmng.FontProperties(fname=fnt)
fp.set_size(9.0)
import mutils.plotting as mplt
colors_ = mplt.colorset_distinct
import matplotlib as mpl
mpl.rcParams['mathtext.default'] = 'regular' # this is the key
In [10]:
# --- only "top" subplot
fig = figure(figsize=(8.8/2.56,2.8)) # height was 2.4
ax0 = fig.add_subplot(1,1,1)
posP = arange(rvpr_full.shape[1]) * 3
b1 = boxplot(rvpr_full, positions = posP + .1, widths=.30, sym='')
b2 = boxplot(rvpr_fac, positions = posP + .5, widths=.30, sym='')
b3 = boxplot(rvpr_ank, positions = posP + .9, widths=.30, sym='')
b4 = boxplot(rvpr_aug, positions = posP + 1.3, widths=.30, sym='')
mi.recolor(b1, colors_[0], lw=.5)
mi.recolor(b2, colors_[1], lw=.5)
mi.recolor(b3, colors_[2], lw=.5)
mi.recolor(b4, colors_[3], lw=.5)
fmt = {'marker' : '>',
'markersize' : 3,
'markeredgecolor' : 'None'}
for idx in nvs:
for x in arange(5):
plot(posP[x] + .1 - .1, median(a_rvpr_full[idx][:,x]), color=colors_[0], **fmt)
plot(posP[x] + .5 - .1, median(a_rvpr_fac[idx][:,x]), color=colors_[1], **fmt)
plot(posP[x] + .9 - .1, median(a_rvpr_ank[idx][:,x]), color=colors_[2], **fmt)
plot(posP[x] + 1.3 - .1, median(a_rvpr_aug[idx][:,x]), color=colors_[3], **fmt)
ax0.set_xticks(arange(rvpr_full.shape[1]) * 3 + .7)
ax0.set_yticks(arange(6) * .2)
ax0.set_yticklabels(arange(6) * .2, fontproperties=fp)
ax0.set_xticklabels(['$k$', r'$\alpha$', r'$L_0$', r'$\beta$', r'$\Delta E$'], fontproperties=fp)
ax0.set_xlabel('predicted parameter', fontproperties=fp)
ax0.set_ylabel('rrv', fontproperties=fp)
ax0.plot([-1, rvpr_full.shape[1] * 3 - 1], [1, 1], 'k--')
ax0.set_xlim(-2, rvpr_full.shape[1] * 3 - 1)
ax0.set_ylim(0, 1.15)
subplots_adjust(right=.966, bottom=.16, top=.77)
ax0.set_title('Predictability of SLIP parameters\n\n\n', fontproperties=fp)
axl = fig.add_axes([0.12, 0.79, 0.89, 0.097], frameon=False )
fp.set_size(8.0)
axl.set_xticks([])
axl.set_yticks([])
axl.text(1, 0, 'full state', fontproperties=fp, va='center')
axl.text(1, -1, 'CoM and factors', fontproperties=fp, va='center')
axl.text(7, 0, 'CoM and ankle state', fontproperties=fp, va='center')
axl.text(7, -1, 'CoM and prev. SLIP parameter', fontproperties=fp, va='center')
axl.plot([0, .8],[0, 0], '-', color=colors_[0], lw=2)
axl.plot([0, .8],[-1,-1], '-', color=colors_[1], lw=2)
axl.plot([6, 6.8],[0, 0], '-', color=colors_[2], lw=2)
axl.plot([6, 6.8],[-1,-1], '-', color=colors_[3], lw=2)
axl.set_ylim(-1.5,.5)
axl.set_xlim(0,15)
fp.set_size(9.0)
ax0.arrow(-1.5,0.95,0,-0.8, head_width=0.15, head_length=0.1, lw=1, )
ax0.text(-1.3,0.5, 'better prediction', ha='left', va='center', rotation=90, fontproperties=fp)
if True:
savefig('img/slip_param_prediction_bxpltS{}.pdf'.format(
vis_subj_str))
pass
In [18]:
# cell to prevent excessive scrolling
In [11]:
close('all')
In [21]:
# --- visualize (only "bottom" subplot)
# --- --- select subject for "full" visualization
vis_subj = 2 # 0,1,2,3
c2_idx = 8 # color index: 2="old", 8="new" (lighter green)
pull_ank_down = True # pull ankle labels "down" on same axis
nvs = list(set(range(4)) - set([vis_subj]))
vis_subj_str = ['1', '2', '3', '7'][vis_subj]
rvpr_full = a_rvpr_full[vis_subj]
rvpr_fac = a_rvpr_fac[vis_subj]
rvpr_ank = a_rvpr_ank[vis_subj]
rvpr_aug = a_rvpr_aug[vis_subj]
rvfr_fac = a_rvfr_fac[vis_subj]
rvfr_full = a_rvfr_full[vis_subj]
rvankr_ank = a_rvankr_ank[vis_subj]
rvankr_full = a_rvankr_full[vis_subj]
# --- --- create figure
if pull_ank_down:
fig = figure(figsize=(18./2.56, 2.8))
else:
fig = figure(figsize=(18./2.56, 3.2))
# --- --- right subplot
outliers = '' # enter '+' to show outliers
ax1 = fig.add_subplot(1,1,1)
pos1 = hstack([arange(3) * 3, arange(conf.n_factors) * 1.8 + 10])
pos2 = hstack([arange(3) * 3, arange(6) * 1.8 + 19]) + .8
b1 = boxplot(rvfr_fac, sym=outliers, positions=pos1 + .1, widths=.40)
b2 = boxplot(rvfr_full, sym=outliers, positions=pos1 + .8, widths=.40)
if not pull_ank_down:
ax1b = ax1.twiny()
else:
ax1b = ax1
b3 = boxplot(rvankr_ank, sym=outliers, positions=pos2 + .9, widths=.40)
b4 = boxplot(rvankr_full[:,3:], sym=outliers, positions=pos2[3:] , widths=.40)
mi.recolor(b1, colors_[1])
mi.recolor(b2,colors_[0],)
mi.recolor(b3, colors_[c2_idx]) # was: colors_[2]
mi.recolor(b4,colors_[0],)
# remove outliers
if pull_ank_down:
xticks1 = hstack([pos1[:3] + .8, pos1[3:] + .35, pos2[3:] + .45])
else:
xticks1 = pos1 + .45
xticks1[:3] += .35
xticks2 = pos2[3:] + .45
ax1b.set_xticks(xticks2)
ax1.set_xticks(xticks1)
ax1.grid('on')
ax1b.grid('on')
#ax1.set_xticklabels(['com$_z$', 'com$_{vx}$', 'com$_{vy}$',
# 'factor 1', 'factor 2', 'factor 3', 'factor 4', 'factor 5'],
# fontproperties=fp, rotation=45, ha='right')
#ax1b.set_xticklabels(['ank$_x$', 'ank$_y$', 'ank$_z$', 'ank$_{vx}$', 'ank$_{vy}$', 'ank$_{vz}$'], fontproperties=fp,
# rotation=45, ha='left')
if pull_ank_down:
ax1.set_xticklabels(['com$_y$', 'com$_{vz}$', 'com$_{vx}$',
'factor 1', 'factor 2', 'factor 3', 'factor 4', 'factor 5',
'ank$_z$', 'ank$_x$', 'ank$_y$', 'ank$_{vz}$', 'ank$_{vx}$', 'ank$_{vy}$'],
fontproperties=fp, rotation=45, ha='right')
else:
ax1.set_xticklabels(['com$_y$', 'com$_{vz}$', 'com$_{vx}$',
'factor 1', 'factor 2', 'factor 3', 'factor 4', 'factor 5'],
fontproperties=fp, rotation=45, ha='right')
ax1b.set_xticklabels(['ank$_z$', 'ank$_x$', 'ank$_y$', 'ank$_{vz}$', 'ank$_{vx}$', 'ank$_{vy}$'], fontproperties=fp,
rotation=45, ha='left')
#ax1b.set_yt
#ax1.set_xlim(-1.5, pos2[-1] + 2)
ax1.set_xlim(-2, pos2[-1] + 2) # additional space for "better prediction" arrow
ax1.arrow(-1.5,0.95,0,-0.8, head_width=0.225, head_length=0.1, lw=1, fc='k', ec='k')
ax1.text(-1.3,0.5, 'better prediction', ha='left', va='center', rotation=90, fontproperties=fp)
ax1b.plot([-2., pos2[-1] + 2], [1, 1], 'k--')
ax1b.set_xlim(-2, pos2[-1] + 2)
ax1.set_yticklabels(arange(6) * .2, fontproperties=fp)
ax1.set_xlabel('predicted coordinate', fontproperties=fp)
ax1.set_ylim(0, 1.15)
ax1.set_ylabel('rrv', fontproperties=fp)
#ax1b.set_title('Consistency check of reduced model: prediction of' +
# ' own states (boxplots for subject {})\n\n\n\n'.format(vis_subj_str),
# fontproperties=fp)
ax1.set_yticks(arange(6) * .2)
for idx in nvs:
for x in arange(len(pos1)):
plot([pos1[x] +.2, pos1[x] +.7],
[median(a_rvfr_fac[idx][:,x]), median(a_rvfr_full[idx][:,x])], '-', color=colors_[1], linewidth=1)
plot([pos1[x] +.2],
[median(a_rvfr_fac[idx][:,x])], '>', color=colors_[1], mec='None',
markersize=5)
plot([pos1[x] +.7],
[median(a_rvfr_full[idx][:,x])], '<', color=colors_[0], mec='None',
markersize=5)
for x in arange(len(pos2)):
plot([pos2[x] + .8, pos2[x] + .1],
[median(a_rvankr_ank[idx][:,x]), median(a_rvankr_full[idx][:,x])], '-', color=colors_[c2_idx], linewidth=.5)
plot([pos2[x] + .8],
[median(a_rvankr_ank[idx][:,x])], '<', color=colors_[c2_idx], mec='None',
markersize=5)
plot([pos2[x] + .1],
[median(a_rvankr_full[idx][:,x])], '>', color=colors_[0], mec='None',
markersize=5)
if pull_ank_down:
subplots_adjust(hspace=.5, right=.975, top=.95, bottom=.24, left=.075)
else:
subplots_adjust(hspace=.5, right=.975, top=.8, bottom=.21, left=.075)
#subplots_adjust(right=.966, bottom=.18, top=.9)
#savefig('tmp/prediction_and_consistency.pdf')
# --- legend axis
if pull_ank_down:
axl = fig.add_axes([0.1, 0.875, 0.75, 0.0625])
else:
axl = fig.add_axes([0.1, 0.78, 0.55, 0.125])
axl.set_xticks([])
axl.set_yticks([])
if pull_ank_down:
axl.text(0, -1, 'input for prediction:', fontproperties=fp, va='center')
xoffset = 15
else:
axl.text(0, 0, 'input for prediction:', fontproperties=fp, va='center')
xoffset = 0
axl.text(3 + xoffset, -1, 'full state', fontproperties=fp, va='center')
axl.text(13 + xoffset, -1, 'CoM and factors', fontproperties=fp, va='center')
axl.text(28 + xoffset, -1, 'CoM and ankle state', fontproperties=fp, va='center')
axl.plot([0 + xoffset, 2 + xoffset],[-1,-1], '-', color=colors_[0], lw=2)
axl.plot([10 + xoffset, 12 + xoffset],[-1,-1], '-', color=colors_[1], lw=2)
axl.plot([25 + xoffset, 27 + xoffset],[-1,-1], '-', color=colors_[c2_idx], lw=2)
if pull_ank_down:
axl.set_ylim(-1.5,-.5)
else:
axl.set_ylim(-1.5,.5)
axl.set_xlim(-2, 42 + xoffset)
savefig('img/consistency_check_prediction_s{}.pdf'.format(vis_subj_str))
print 'stored as img/consistency_check_prediction_s{}.pdf'.format(vis_subj_str)
pass #suppress output
In [17]:
[(nr, val) for nr, val in enumerate(colors_)]
Out[17]:
In [ ]: